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Abstract 

For a particular experimental design, there is interest in finding 
which polynomial models can be identified in the usual regression set 
up. The algebraic methods based on Grobner bases provide a sys- 
tematic way of doing this. The algebraic method does not in general 
produce all estimable models but it can be shown that it yields mod- 
els which have minimal average degree in a well-defined sense and in 
both a weighted and unweighted version. This provides an alterna- 
tive measure to that based on "aberration" and moreover is applicable 
to any experimental design. A simple algorithm is given and bounds 
are derived for the criteria, which may be used to give asymptotic 
Nyquist-like estimability rates as model and sample sizes increase. 



1 Introduction 

It is of considerable value to represent an experimental design as the solution 
of a set of polynomial equations. In the terminology of algebraic geometry a 
design is a zero dimensional variety and the corresponding ideal comprising 
all polynomials which are zero on every design point is called an "ideal of 
points". Pistone & Wynn (1996) first used explicit methods from algebraic 
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geometry and in particular introduced Grobner bases into designs. Issues to 
do with identifiability of polynomial regression models, or inte rpolators, can 



be tra nslated into problems about such varieties and ideals, see iPistone et al 



(2001) 



The purpose of this paper is to introduce the notion of linear aberration of 
a polynomial model. Linear aberration is defined only for polynomial models, 
which are used routinely in statistical literature. A polynomial model with 
low order terms has low aberration, thus engaging low aberration with the 
standard practice of preferring polynomial models with low order terms. The 
preference for mode l s with low order terms has been acknowl edged in recent 



papers, see lLi et al.l (120031 ) and lBalakrishnan and Yangi (120061 ). although they 



do not refer to linear aberration. 

Let a = (cui, . . . , ad) be a nonnegative ci-dimensional integer multi- index. 
A monomial in the indeterminates x\,...,Xd is the power product x a = 
A model basis is a collection of distinct monomials {x a ,a £ L}, 
where L is a finite set of multi-indices. By combining linearly monomials in 
L we form polynomials: 

Vl{x) = y^JgX 01 , 

a<=L 

where 9 a are real coefficients. The polynomial t]l(x) is a candidate for inter- 
polation or statistical modelling. 

This paper is concerned with the following concepts. 

Definition 1 Let L be a model basis and let w = (wi, . . . , wj) be a collection 
of non-negative weights with ^ =1 «jj = 1. We define the weighted linear 
aberration of L as 



Wja 



n 

(ai,...,a d )GL i=l 



where n is the number of elements in L. 

We are interested in studying aberration for models identifiable by an 
experimental design and along this paper we compare models and designs of 
the same size n. 

Definition 2 An experimental design D, of sample size n = \D\, is a set of 
points in IR d . 
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We say that a model basis L with cardinality \L\ = n is identifiable by D 
if the design model matrix X = [x a } x& E>,a&L is invertible. 

The term aberration is used to acknowledge the work on "minimum 
aberration" for r e gular fracti onal factoria l desig ns of Wu and others, see 
Fries and Hunter (jl980h and IWu and Wul (120021). For fractional factorial 



designs, the notion of estimation capacity is related to the abil i ty of a de- 
sign to identify models of low degree, see ICheng and Mukerjed (119981 ) and 



Chen and Cheng! (120041 ) . We do not make a direct mathematical comparison 



with that work but simply point to a common motivation. 

In Section [2] we review the basic ideas on algebraic identifiability. The 
search for identifiable models is driven by a divisibility condition, which 
makes the search problem tractable. We then introduce the state polytope, 
whose vertices correspond to the models identified using the algebra. In 
Section [3] we study aberration. The basic ideas on aberration are closely 
linked with the algebraic work on corner cut models and state polytopes in 



Onn and Sturmfelsl (119991 ) . We are specially interested in obtaining minimal 
values for aberration for which we establish upper and lower bounds. An 
approximate approach to minimal aberration is discussed. In Section H] we 
discuss various examples. In Section [5] we discuss possible extensions of the 
theory and, by example, a connection with the notion of aberration by Wu 
and others is discussed. 



2 The G-basis method and the state polytope 



The aberation A(w, L) has remarkable connections with t he algebraic method 
i n exp erimental design introduced by Pis tone and Wynn iPistone and Wynn 
( 1996 ) and developed in the monograph Pistone et al.l ( 2001 ) and the joint 
work of Onn and Sturmfels Onn and Sturmfelsl (Il999f ). In this section we 
present the basic ideas on identifiability using algebraic techniques. 

Let the set of all monomials in d indeterminates be T d = {x a , a G Z> }, 
where Z> is the set of non-negative integers and Z> is the set of all vectors 
in d dimensions and with entries in Z> . A polynomial is a finite linear com- 
bination of monomials in T d with real coefficients. The set of all polynomials 
is denoted as M[xi, . . . ,Xd\. It has the structure of a ring with the usual 
operations of sum and product of polynomials. 

A term ordering y on R[xi, . . . , xa] is a total ordering on T d such that i) 



x 



y 1 for all x^r.a/ffl,..., 0) and ii) for all x c 



e T d if x a y x 13 
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then x a x~ 1 y x^x 1 . The leading term of a polynomial is the largest term with 
non-zero coefficient with respect to y. For a polynomial / G R[xi, . . . ,xj\, 
we write its leading term as LT^(/). 



A partial order on T is defined by a vector w G 



as x a >~ w x 13 if 



a > w T j3, where x a , x 1 3 G T d and w T is the t ransposed vector of w. Under 



Babson et all (l2003h ; ICox et all (Il997f ) ) this defines 



w 

some conditions on w [see 
a term order. Given a term order there are w such that x a y x 13 if and 
only if x a y w x 13 . 

A design D, considered as a zero-dimensional variety gives rise to a de- 
sign ideal, 1(D), which is the set of all polynomials which have zeros at 
all the points of D. We have that 1(D) C M[xi, . . . , x^]. The polyno- 
mial ideal I is generated by the set of polynomials G = {gi,...,g s } if 
1 = {Yh=i fi9i '■ fi e ^[ x u ■ ■ ■ , x d]} and we write / = (g x , ...,g s ). 

An important set of generators for the desig n ideal is the Grob ner basis. 
Grobner bases were introduced by Buchberger in lBuchbergerl (|19661) and they 
have become a powerful computational tool in many fields I Cox et al.l (119971 . 
20051 ). A Grobner basis of 1(D) with respect to a term order y is a finite 



subset G y (D) c 1(D) such that (LT y (g) : g G G y (D)) = (LT^(f) : f G 
1(D)). The computation of Grobner bases is implemen ted in standa r d com - 
puter programs suc h as CoGA, Singular or Maple, see ICoCoATeaml (120071 ); 
Greuel et all fl2005f ): iMonagan et all (j2005h . 

Two polynomials / and g in . . . , Xj] are equivalent with respect to 
1(D) if the following conditions hold: 

i) f-gel(D) 

ii) f(d) = g(d) for all d G D 

Given a term ordering >-, the quotient ring M[xi, . . . , Xd\/I(D) has a unique 
M-vector space basis given by the monomials in T d that cannot be divided 
by the leading terms of the polynomials in G y (D) for 1(D). The monomial 
basis so obtained, or equivalently, the set of its exponents L = L(D,y), 
has a staircase (also echelon, order ideal) property: for a G L, if (3 < a 
componentwise, then f3 G L. Equivalently we say that for any x a G L, if 
x 13 divides x a then x 13 G L. We call bases which have a staircase structure 
staircase models. The dimension of WL[xi, ... ,Xd]/ 1(D) as R-vector space 



is n, see iPistone and Wynnl (119961 ). i.e. the number of points in D and of 



multi-indices in L is n. 
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For a given basis of the quotient ring with exponents in L and a set of 
real values (data) Y x ,x G D, there exists a unique interpolator r] L (x) such 
that = T]l(x), x G D. Other non-saturated statistical sub-models can be 



constructed from subsets of L, see iHollidav et ail fll999h and lPeixotol fll987h 



Definition 3 The algebraic fan of D is C a (D) = {L(D,y), where y is a 
term ordering in M[xi, . . . ,Xd}}- This is the collection of staircases L(D, y) 
arising from a fixed design D by varying all monomial orderings. 



The algebraic fan of a design was proposed by Caboara et al. ICaboara et al. 



(199 7|), constructing upon t he algebraic fan of an ideal of Mora an d Robbiano 
Mora and Robbiand (119881 ) . Babson et al. iBabson et al.l (120031 ) proposed a 
polynomial time algorithm to compute C a (D). They compute an efficient 
set of weight vectors and perfo rm a change o f basis which stems from the 
so-called FGLM algorithm, see iFaugere et al.l (119931 ). In Section [3.11 an al- 
gorithm is presented to identify a model in the algebraic fan using a weight 
vector. 

It is important to note that not all staircase models identified by D are 
in C a (D). We denote the set of all identifiable staircase models for a design 
D as C S (D). In fact the algebraic fan is small relative to C S (D), that is 
£>a(D) C C S (D), see Chapter 6 in the unpubl i shed P h.D. thesis by Maruri- 
Aguilar (2007) and Section 4 in lPistone etaD (I2006T ). 

We now establish the link between the algebraic fan of a design and the 
state polytope of the design ideal. For a model basis L = {a 1; . . . , a n }, G 



Z> define 



OtL 



E 



0£i. 



This vector appears in the definition of A(w, L) and we can write A(w, L) 



:w 



T a>L)/n. The set all such vectors over C a (D) gives the state polytope. 



Definition 4 The state polytope S(D) of a design D, or equivalently of the 
design ideal 1(D) is the convex hull 

S(D) := conv ({«£, : L is a staircase in C a (D)}) . 



The following theorem (jSturmfeld . 1 19961 . Ch. 2) summarizes the connec- 
tion between the state polytope and the set of models C a (D), i.e. the relation 
between a design and its algebraic fan. 
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Theorem 1 Let D be a design and let S(D) be its state polytope. Then the 
set of vertices of the state polytope of D is in one to one correspondence with 
the algebraic fan of D. 

The state polytope does not only contain information concerning models 
in the algebraic fan of a design, but it also provides information about the 
term ordering vectors needed to construct it. We recall that a (^-dimensional 
polytope is a bounded subset of M. d , which corresponds to the solutions of a 
system of linear inequalities. The normal cone of a face of a polytope is the 
relatively open cone of those vectors in M. d uniquely minimised over the face 
of the polytope. The normal fan of a polytope is the collection of all the 
normal cones of the polytope. 

Two ordering vectors w and w' are said to be equivalent (modulo 1(D)) 
if L(D, y w ) = L(D, y w >). The normal fan of the sta te polytope pa r tition s 



into equiyalenc e classes of orderi ng vectors, see iBabson et al.l (120031 ); 



Fukuda et ah] (120071 ); ISturmfelsl (119961 ) . Indeed every vertex of S(D) corre- 



sponds to a model in £ a (D). Moreover, the interior of the normal cone of 
a vertex in S(D) contains those vectors w which correspond to the same 
equivalence class. 

We motivate Theorem [2] below with a simple example. The black dots 
in Figure [1] give a 5 point design in 2 dimensions, D. They also give the 
set of exponents L obtained for any term ordering, indeed the size of the 
algebraic fan of D is one. The crosses represent the exponents of the leading 
terms of the Grobner basis: (2, 0), (1, 2), (0, 3). The line separates the model 
exponents, L, from these leading terms. This is an example of a corner cut 
model. Note that equivalently the line separates L from its complement in 

Definition 5 A model L, of size \L\ = n, is said to be a corner cut model if 
there is a (d — 1) dimensional hyperplane separating L from its complement 
Z|o\L. 

Not all staircases are corner cuts, for example L = {(0, 0), (1, 0), (0, 1), (1, 1)} 
is a staircase that cannot be separated by a hyperplane from its complement 
in Z| . 

The set of exponents of a corner cut model is referred to as a corner cut 
staircase or simply, as a corner cut. Corner cuts were introduced by Onn and 



Sturmfels lOnn and Sturmfelsl (119991 ). A generating function for the number 
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Figure 1: Corner cut and separating hyperplane. 



of bidimensional corner cuts is given in ICorteel et al.l (119991 ). while the order 
of the card i nality of the set of corner cuts is proven bounded by {n\ogn) d ~ l 



m 



Wagnerl (120021 ) . A special class of designs is composed with those designs 
that identify all corner cut models of a given size. 

Definition 6 A design D C K d comprised of n distinct points is said to be 
generic if all corner cut models of size n = \D\ are identifiable. 

A special polytope is constructed with the exponents for corner cut mod- 
els. It will be used to compute the algebraic fan of generic designs. 

Definition 7 The corner cut polytope is CC(n,d) := conv({c*j, : L is a 
corner cut staircase in d dimensions and of size n}). 

For a discussion on th e prope r ties o f bidimensional corner cut polytopes 



see the paper by Miiller iMuilerl (120031 ). The algebraic fan of generic de- 



signs corresponds to the set of corner cut models, as stated in the following 
theorem. 

Theorem 2 (Onn and Sturmfels, 1999) Let D C M. d be a generic design 
with n points. Then 

i) S{D) = CC{n,d) and 

ii) the algebraic fan of D is the set of corner cut models in d dimensions 
and with n elements. 

We remark that the corner cut polytope is an invariant object for the 
class of all the ideals generated by generic designs with the same sample 
size n and number of factors d and all generic designs have the same state 
polytope. 
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3 Minimal linear aberration 



An important feature of the state polytope is that its vertices are automat- 
ically "lower" vertices in the sense of convexity. State polytopes relate di- 
rectly to models with minimal linear aberration. In Section I3TT1 an algorithm 
to compute a models of minimal aberration is presented. 

Theorem 3 Given a design D C M. d with n distinct points and a weight 
vector w G 1R> , there is a least one vertex a* G S(D) which minimises 
A(w,L) over all identifiable staircase models C S (D), that is 

-(w T a*) = A(w,L*) = min Aiw.L) 
n LeCs(D) 

for all L* such that ctL* = ct* . Moreover, given a vertex of S(D), there is at 
least one w* G M> such that this vertex (model) minimizes A(w,L), that is, 

A(w*,L) = min A(w, L) 



for L such that ot^ = a^. 

Proof. First, for given w we minim ise w T aL for L G C a (D), which is a finite 



set, see 



Mora and Robbiand (jl988l ). The ql for L G C a {D) are vertices of 



S{D) by definition. Furthermore, because we rest rict L to the alge braic fan 



of D there cannot be three aligned o7^ in S(D), see ISturmfelsl (119961 ). For the 
second claim, it is sufficient to take a vector wl in the interior of a normal 
cone for a^. By definition, A(w,L) is minimised for vectors on the interior 
of the normal cone. ■ 

Theorem H] follows directly from Theorem [31 

Theorem 4 For every weight vector w there is a design D C M. d which 
minimizes A(w,L), among all designs with sample size n and identifiable 
staircases. 

This is stated compactly as: 

A*(w,n) = min min A(w,L) 

D:\D\=n L&C a (D) 

is achieved for a generic design. That is, if a design is generic then automat- 
ically its algebraic fan contains models of minimal aberration. 
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3.1 Computation of the minimal aberration model 

The model minimizing linear aberration can be found by the greedy algorithm. 
Let D be a design; let w be a fixed weight vector in M> and let T be the 
following set of potential exponents 

r := la = (ai, . . . , a d ) G Z|„ : + 1) < n 



The set T contains all staircase models with n terms, see lBabson et al.l (120031 ). 
Now define the weight of a G rtobea;(a) := -J2i=i w i a i = (w T a)/n. Order 
the vectors in T by their weight uj(-) in increasing order, that is, index them 
as a 1 , . . . , ct' r ' such that ^(a 1 ) < • • • < cu(a^), where |T| is the cardinality 
of T. Then the set L C T with the first n terms of T which are identifiable 
by D has minimum aberration. 

The model basis L is constructed by the following procedure: initialize 
L := 0; while \L\ < n, find a % of smallest index with respect to u(-) such that 
the column vectors d a , a G LU{a 1 }, d G D, are linearly independent; update 
L := L U {a 1 } and repeat until \L\ = n. We have the following theorem. 

Theorem 5 Let D C M. d be a design; let w be a fixed weight vector with pos- 
itive entries and let L be the model basis constructed by the greedy algorithm. 
Then L belongs to the algebraic fan of the design. 

Example 1 Consider the design D = {(0,0), (1,0), (0, 1), (-1, 1)} and the 
weight vector w = (4, 1). The set of potential exponents, T contains 8 ele- 
ments, which are sorted out using the weight function u(-) as 

T ={ (0,0), (0,1), (0,2), (0,3), (1,0), (1,1), (2,0), (3,0) } 
nuj(-) = 0123458 12 

The first 4 elements in T such that their design columns are linearly indepen- 
dent are L = {(0,0), (0, 1), (1,0), (0, 1)}. Thus the set L of minimal linear 
aberration corresponds to the model with terms {1, Xi,X2, £1X2}. 

3.2 Examples 

We can compare different designs using aberration as long as they have the 
same number of factors d and the number of points n. For a design D, the 
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state polyhedron of D is obtained by ( Mink owski) addition of M> to the 



state polytope S(D), see iBabson et al.l (120031 ). The state polyhedron yields 
the same information as the state polytope. Indeed the normal fan of the 
(nega tive) state polyhedro n yields automatically the first orthant, see Fukuda 



et al. iFukuda et all (120071 ). 
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Figure 2: The left graph depicts S(D) and the state polyhedron for the CCD 
of Example [2j The right graph shows state polyhedra for the three designs of 
Example [21 The empty dots correspond to vertexes/models identified by the 
generic design only, while the triangle is for the sole model in the algebraic 
fan of the 3 2 design. 



Example 2 Co n sider a central composite design (CCD by Box and Wilson 
Box and Wilson! (jl95ll )) with two factors, one observation at the origin and 
axial distance a = a/2. The CCD has 9 runs and its algebraic fan contains 
exactly two models, namely 



{1, Xi, xl, x\, x\, 2 2 , XxX2, x\x 2 , xQ 



(1) 



together with the model obtained by permuting the roles of x\ and x^- Let 
L\ be the set of exponents of the model support in Equation (pQ). Clearly, 
a^i = (13, 5) and the state polytope for the design ideal of the CCD is 
conv ({(13, 5), (5, 13)}), see left graph of Figure El No w consider a g eneric 



desi gn with the same number of runs as the CCD. In ICorteel et al.l (119991 ) 



and lOnn and Sturmfeld (119991 ) it is shown that there are 12 corner cut models 



for d = 2 and n — 9. By Theorem [21 the algebraic fan of the generic design 
contains all the 12 corner cut models, including those in the algebraic fan 
of the CCD. We consider also a full factorial design 3 2 , which identifies only 
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the model with support {\,X\,x\} <8> {l,X2,x%}, where <8> is the Kronecker 
product. Its state polytope is the point (9, 9). In the right graph of Figure [2] 
we depict the state polyhedra for the three designs and in Figure [3] we plot 
mm Le£ a (-D) A(w, L) for w = (w\, w^) € [0, l] 2 and W\ + W2 — 1. For the CCD 
this is 

f ((wi,l-wi)(13,5) T )/9 = (8wi + 5)/9 if w x < 1/2 
\ ((wi,l -Wi)(5,13) T ) /9 = (-8101 + 13)/9 if w x >\/2 

For the generic design the aberration curve is a piecewise linear function 
with 12 segments. Finally, the aberration for the design 3 2 is constant for all 
weights. As expected, the aberration takes its minimum value for the generic 
design, over all possible weights. 




Figure 3: Minimal aberration for three designs in two factors and nine runs, 
see Example [21 



Example 3 Consider the design D = {(0,0), (1, 1), (2, 2), (3, 4), (5, 7), (11, 13), 
(a, (3)}, where (a, ft) « (1.82997,1 .82448) is the only real so lution of a sys- 
tem of polynomial equations, see (jOnn and Sturmfeld . fl99~9l . Page 47). The 
algebraic fan of the above design has ten models and its state polytope is 

conv ({(21, 0), (15, 1), (11, 2), (9, 3), (6, 5), (5, 6), (3, 9), (2, 11), (1, 15), (0, 21)}) . 

Now consider a generic design G with the same number of runs and factors. 
The algebraic fan of G is the set of corner cut models which for 7 points in 2 
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dimen sions has 8 elements, see ICorteel et al.l (119991 ) and lOnn and Sturmfels 
(119991 ) and thus its state polytope is the corner cut polytope: 



CC(7, 2) = conv ({(21, 0), (15, 1), (11, 2), (7, 4), (4, 7), (2, 11), (1, 15), (0, 21)}) . 

In Figure H] we graph the aberration for both designs as a function of w\. 
Although the size of the algebraic fan of D is bigger than that for a generic 
design, the weighted aberration takes minimal value for the generic design 
for all possible weight vectors (w\, 1 — w\). 




Figure 4: Minimal aberration for G (solid line) and D (dashed line), see 
Example [3l 



Example 4 The aberration of some sets of multi-indices does not depend 
on w. For instance, consider the following sets in two dimensions 

L n = : i = 0, . . .,n - 1} 

M n = {(i,j):i,j = 0,...,n-l} 
N n = {(i,j):0<i + j<n} 

for which the aberrations are A(w,L n ) = (n — l)/2, A(w,M n ) = (n — l)/2 
and A(w,N n ) = n/3. To properly compare the above aberrations, the sets 
L, M, N must have the same size. Below are values m, n and p = m 2 such 
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that #L P = #M m = #N n for m up to 8000. 



m 


n 


A(w,L p ) 


A(w, M m ) 


A(w,N n ) 


1 














6 


8 


17.5 


2.5 


2.6 


35 


49 


612.0 


17.0 


16.3 


204 


288 


20807.5 


101.5 


96.0 


1189 


1981 


7.0 x 10 5 


594.0 


560.3 


6930 


9800 


2.4 x 10 7 


3464.5 


3266.6 


40391 


57121 


8.1 x 10 s 


20195.0 


19040.3 



As sample size grows, the aberration of the triangular set N n remains smaller 
than for the square set M m . 

3.3 Bounds for the aberration 

Although the minimal value of the aberration A*(w,n), depends on the 
weight vector w = (wi, . . . ,Wd), we can carry out a special normalisation 
which leads to bounds for the minimal aberration. These bounds depend 
only on a simple function of the weights, surprisingly the geometric mean. 
Our construction is based upon the expected value of auxiliary random vari- 
ables which are suitably constructed. 

For the rest of this Section let D C M. d be a generic design with n points. 
Let w be a fixed weight vector with positive elements and let L be the corner 
cut model identified by w. We recall that \L\ = n. 

For an integer multindex a define its upper cell as the unit cube with 
lower vertex at a 

c{a) = {v G R d : on < v t < an + 1} 
and similarly the lower cell of a is 

c(a) = {v E R d : ai - 1 < Vi < ai] 

Define: 

<2 = U aGL c(a), Q = U Q6L c(a). 

See Figure [5] for a depiction of lower and upper cells with L a corner cut. 

Clearly, the volume of Q and of Q equals n, that is the cardinality of L. 
We now create a simplex S(w) C M d which is directed by the vector w and 
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Figure 5: Bidimensional corner cut together with upper (left diagram) and 
lower cells (right diagram) Q and Q. In both diagrams the vector w, a 
separating hyperplane and equivalent simplexes S(w) and S_(w) were added. 



has volume n. We call this simplex and the subset of the first orthant below it 
the equivalent simplex, which is formally S(w) = |t> £ M> : Y2t=i 
The volume of S(w) is determined up to the constant c > 0. We find the 
value of this constant by setting the total volume of the equivalent simplex 
equal to n: 

c d 

n = ^ j 

d\ n i=1 Wi 

giving 

c={nd\y*g{w), (2) 

where 

/ d 

g(w) = I Y[wi 
\i=i 

is the geometric mean of the components of the weight vector w. We call 
H(w) the hyperplane which limits the equivalent simplex, that is H(w) = 

|u £ Ri : Y? i= i v i w i = c }- 

The expected value of a random variable with uniform support over S(w) 
will be used now to compute bounds for aberration. We can compute a no- 
tional value of A, the linear aberration for a distribution D as the expectation 
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A(w, S(w)) = E(J^ WiXi) for the random vector (Xi, . . . , X d ) with uniform 
distribution over S(w). Thus for the equivalent simplex we have that 



after substituting Equation (J2J) in A(w, S(w)) 

We observe that the region Q is obtained from Q by a negative shift 
(— 1, . . . , — 1). As before, we consider a random vector with joint uniform 
distribution over Q. We then use the expected value of WiXi as the aber- 
ration A(w, Q). Analogously we define A(w, Q) and we have 

A{w,Q)=A{w,Q)-l 

Similarly we can create a region S_(w) by the same downward shift, and we 
have 

A(w,S(w)) = A(w,S(w)) - 1. 

As D is generic and thus L is a corner cut there exist cutting hyperplanes 
separating L from its complement in Z> . Moreover if w is in the interior 
of the normal cone of the corner cut polytope, then we can select a cut- 
ting hyperplane H which is orthogonal to w and thus parallel to H(w), see 



Onn and Sturmfeld (119991 ) . 



Example 5 Consider a generic design with d = 2,n = 3 and L = {(0,0), 
(1, 0),(2, 0)}. The weight vector w = (1, 2) is not in the interior of a normal 
cone of the corner cut polytope CC(2, 3). Indeed the weight vector is on 
the boundary of the normal cone separating L from the corner cut model 
{(0, 0), (1, 0), (0, 1)}. The hyperplanes perpendicular to w are 2x\ — %2 = c 
and none of them is a cutting hyperplane for L. 

By a simple argument the simplex Sh with faces x% — 0, (i — 1, . . . , d) 
and H lies wholly within the upper quadrant region Q because otherwise, 
the cutting hyperplane hypothesis for H would be violated and thus Sh has 
volume less than n. Recall that the equivalent simplex S(w) has volume n. 

There is one additional argument that leads to our first inequality. Since 
the region Q and the equivalent simplex S(w) have the same volume n, it 
must be that Q protrudes beyond S(w). Equivalently we may move mass 
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from Q, that is, beyond H(w), inside S(w). As this mass occurs orthogonally 
to w, we claim that this movement diminishes the aberration, thus 

A(w,S(w)) < A(w,Q). 

This property is also inherited by the downward shifted version, and we have 
A(w,S(w)) < A(w,Q). The same orthogonality argument shows the middle 
inequality in the following sequence: 

A(w,S(w)) < A(w,Q) < A(w,S(w)) < A(w,Q). 

By Theorem HI as the design is generic and L is the model identified by w, 
clearly we have 

A(w,Q) < A*(w,n) < A(w,Q). 

Analogous argument and construction as above shows that A(w, Q) < A(w, S(w 
1. 

Theorem 6 Let D C M. d be a generic design with n points; let w G M. d be 
a vector of positive weights. Then the minimal aberration A*(w,n) satisfies 
the bounds 

A(w, S(w)) - 1 < A*(w, n) < A{w, S(w)) + 1, (4) 

where A(w, S(w)) is computed in Equation (T^j. 

There are various kinds of asymptotic that this formula leads to. From 
the inequality between geometric and arithmetic mean we have g{w) < 4. 
This suggests the condition lim,i MOO g{w) = c/d for some constant < c < 1. 
Now for Wi = (1 + Si)/d, with 5i = 0, and assuming convergence of 5f 
and n = k d , we use use Stirling's approximation to obtain 

lim A*(w, n) = — . 

Such limits may be considered as asymptotic identifiability rates, analogous 
to the more familiar Nyquist rates in Fourier analysis. 

Example 6 For small d and n the bounds of Equation (j3J) are rather coarse. 
Figure [6] shows the bounds A(w, S(w)) ± 1 of Theorem [6] together with the 
minimal aberration A*(w, n), plotted as function of for d = 2 and n = 4. 
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Notice that, as function of w, the minimal aberration A*(w,n) is a piece- 
wise linear graph (this is a general fact, consequence of Definition [T]) , each 
segment corresponding to a different vertex (different corner cut) of the corner 
cut polytope. Figures [7J and [8] give the bounds and minimal aberration for 
n = 20 and n = 100. In Figures O [7] and [8] we also added a curve for the 
approximate aberration which is presented in Theorem [7J below. 




Figure 6: Minimal aberration A*(w,n) (solid line) for a generic design with 
d — 2, n — 4; bounds A(w, S(w)) and A(w, S(w)) ± 1 of Theorem [6] (dashed 
lines). We also show approximate aberration A using Theorem [7J (thin dashed 
line) . 



d 



3.4 Approximated state polytope for generic designs 

Note that as w changes the hyperplanes H(w) are tangent to the surface 
defined by 

n* = ' = «»{if 

i=i v 7 

and the (normalised) centroids of the equivalent simplices lie on the surface 
defined by 

i=l v / 
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Figure 7: Minimal aberration A*(w,n) (solid line) for a generic design with 
d = 2, n = 20; bounds A(w, S(w)) and A(w, S(w)) ± 1 and (dashed lines) of 
Theorem [61 The figure also shows approximate aberration A of Theorem [7] 
(thin dashed line) which almost overlaps the solid line. 



We can solve an equivalent optimisation problem to the computations of 
A(w, S(w)) in terms of the tangent surfaces: for all centroids lying above or 
on the surface of Equation (|3j), the minimum value of A(w, S(w)) is achieved 
at the centroid of the tangent. 

In the above argument, we are essentially using the surface in Equation (jSJ) 
to approximate the lower border of the state polytope for a generic design, 
i.e. the lower border of the corner cut polytope. In order to improve the 
bounds given in Theorem [61 it seems natural simply to take a surface defined 
by 

d 

J]>* + a) = b (6) 

t=i 

with fixed a, b. In Theorem El we have a = ±1 and b = b + in Equation (j5J). 
In Appendix B we discuss an approach to select the values a, b to obtain a 
good approximation of the corner cut polytope. 

The following theorem estimates minimal aberration for generic designs 
using the approximation of Equation (jSJ). The proof is based on simple ideas 
of constrained optimization, see Appendix A. 
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Figure 8: Minimal aberration A*(w,n) (solid line) for a generic design with 
d = 2, n = 100; bounds A(w, S(w)) and A(w, S(s)) ± 1 (dashed lines). The 
approximate aberration A of Equation (J7J (thin dashed line) is also plotted, 
but is undistinguishable from the minimal aberration. 



Theorem 7 Let w = (wi, . . . , Wd) be a fixed positive weight vector; let D C 
M. d be a generic design with n points. Let the state polytope of 1(D) be 
approximated by Equation Then the value 

d 

A(w) = db 1/d g(w) - a ^ w l (7) 

i=i 

is an approximation of A*(w,n). 

We recall that g(w) is the geometrical mean of the components in w. Figures 
El [7] and M give examples (d = 2 factors, n = 4, 20, 100) of the minimal 
aberration A(w) in Theorem [71 The values a, b for each case were selected 
using the technique in Appendix B. 

4 Examples 

In this section we discuss through extended examples other possible uses 
of the ideas on generic designs and aberration. In Section 14.11 we explore 
and conjecture the existence of generic designs over Latin hypercubes for 
all factors and sample sizes. In Section 14.21 we compare fractional factorial 
designs through their state polytopes. 
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4.1 Latin hypercube design 



Latin hypercube designs (LH) were first proposed by McKay et al. iMcKay et al 



(119791 ) in the context of computer experiments. Latin hypercubes are designs 
with reasonable space filling properties and good projections in lower dimen- 
sions. 

Theorem H] relates minimal aberration to generic designs, i.e. if the design 
is generic, then it identifies models of lower weighted degree (and minimal 
aberration) for any weight vector w. In what follows we study LH using 
Definition [6] of generic designs. 

The construction of a Latin hypercube design can be summarised as fol- 
lows. 

1. Divide the range of each factor into n equal segments. 

2. Select a value in each segment using a random uniform distribution, or 
any other continuous distribution. 

3. Randomly permute the list for each factor. 



By Theorem 30 in iPistone et al.l (120011 ). a Latin hypercube design con 



structed as above is generic with probability one. 

We now consider a special type of LH designs. This type is constructed 
by selecting a fixed value in every segment in Step [2j For instance, we could 
select the minimum, maximum or the midpoint value for every segment. 

There are a few obvious cases of LH designs which are not generic, for 
example when the points of the design lie on a line. We have performed 
exhaustive search for a few cases of LH in two dimensions. Our search points 
out to the existence of generic LH for different values of d, n. In fact for 
the values we tried the proportion of generic LH tends clearly to one. See 
Figures M and [TD] for a depiction of the results, where we additionally plot the 
proportion of maximal fan designs among LH, i.e. LH designs that identify 
all possible staircase models for given d, n. We have the following conjecture 
for the existence of generic LHS for any value of d, n. 

Conjecture 8 For every d > 2 and n > 2 there exists at least one generic 
LH design, constructed by setting a fixed value for every one of the n segments 
in the above procedure. 
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Figure 9: Percentage of generic LHS designs for d = 2 and n < 15. 
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Figure 10: Minus logarithm of the percentage of non generic LHS designs 
for d = 2 and n < 15. 
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Figure 11: LH on [0, l] 2 for d = 2, n = 10 which are not generic and identify 
L. 



Experimentally we observed that when the sample size is n = f^ 1 ) for 
k > 1, the genericity of a LH design is closely linked to the identification of 
a model of total degree k — 1. For example for k = 4,d = 2,n = 10 there are 
10! LH of which 99% are generic. Of the remaining 1% which are not generic 
only 6 designs (up to reflection and rotation), which are given in Figure [TT| 
identify the cubic model with exponent set 

L = {(0, 0), (1, 0), (0, 1), (2, 0), (1, 1), (0, 2), (3, 0), (2, 1), (1, 2), (0, 3)}. 
4.2 Orthogonal fractions 

In this Section we consider some of the techniques of this paper for the 
class of fractional factorial designs with two levels. We first explore the 
relation between state polyhedron and then later propose a tool to compare 
the identification capability of designs. 

In Examples [2] and [3] we observed that in general, nesting of state poly- 
hedra for two designs does not imply any easy relation between the algebraic 
fan of the designs. If instead we restrict to the family of designs with two 
levels then there is a clear relation between such nesting and algebraic fans. 
We have the following Lemma from Chapter 6 in the Ph.D. thesis by Maruri- 
Aguilar (2007). 
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Lemma 9 Let F\ and F 2 be two fractional factorial designs with two lev- 
els and let S\ and S 2 be their corresponding state polyhedra of I (Fi) , I (F 2 ) . 
Then the nesting of state polyhedra S\ C S 2 implies nesting of algebraic fans 
£.(fi)c£,(F a ). 

The following example is based upon Lemma O and presents an interesting 
relation between resolution and identifiability. That is, bigger resolution 
points to more models in the algebraic fan. 

Example 7 Let F\ and F 2 be the 2^ and fractional fractional designs 
with eight runs in four factors and respective generators x\X 2 x^x^ — 1 = 
and X\X 2 Xs — 1 = 0. The subindic es III, IV refer to the resolution of the 
fraction, see Box and Hunter fll961al lbh. Their corresponding state polyhedra 



are nested, i.e. S(F 2 ) C S(Fi) and by direct computation we confirm that 
the algebraic fans are also nested. The algebraic fan £ a (F 2 ) has four models, 
while C a (Fi) includes 12 elements. 

For fractional factorial designs, the estimation of interactions in a design 
was related to the re solution of the design through the property termed 



hidden projection, see lEvangelaras and Koukouvinosl (120061 ) ; IWang and Wu 



( 119951 ). We conjecture the nesting of algebraic fans of two designs 2 k ~ p with 
different resolution. However, exploiting this nesting property of fans to 
compare designs using aberration might need additional considerations. 

Example 8 Let F\, F 2 be the fractions 2 IV ~ 2 given by generators xq— xix 2 x^ = 
0, xj — x 2 xsX4 = and xq —x\x 2 x?,Xi = 0,xi — x\x 2 x^,x^ = respectively. Al- 
though both fractions have the same resolution, the fracti on F 2 corresponds 



to a minimum aberration design using the definition of iFries and Hunter 



([19801) . The state polyhedron 5(F a ) has 133 vertices while S{F 2 ) has 1708. 



There is no nesting of the state polyhedra and C a (Fi) n C a (F 2 ) ^ 

A proposal to compare two designs D\ , D 2 of the same size through their 
state polytopes is to map the vertices of the state polytopes S(Di), S(D 2 ) 
with a function / : M. d — > M.. In this way the state polytopes of D\ and D 2 
are compared by the univariate projections of their vertices. We propose a 
weighted sum of the vertex coordinates 

d 

f(vi,...,v d ) = ^WiVi, (8) 
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with positive weights Wi > 0. We use lOj = 1 for i = 1, . . . , d and thus 
Equation (jSJ) allows for direct comparison of designs based on the distribution 
of total degrees for models in the algebraic fan. 

Example 9 (Continuation of Example [8]) We transform the vertices of the 
state polytopes for F\ and F 2 using Equation (jSJ). In Tabled] in Appendix B 
we summarize the results for each fraction as the distribution of absolute and 
relative frequencies. Clearly, the fraction F 2 with minimum aberration for 
generators identifies models with a smaller total degree than that for F\ and 
in that sense it has smaller linear aberration. See Figure [12] for a histogram 
of the relative frequencies for F± and F 2 . 



-F, 



Figure 12: Histograms of relative frequencies for fractions F\ and F 2 , see 
Example [9] We added F 3 of Example [TO] 



5 Discussion 

5.1 Generalised concave aberration 

This paper is partly concerned with a problem of linear programming, i.e. 
optimising a linear function / : W d —>■ R over a convex polytope. We now 
discuss extensions of our work using other types of aberration. When we 
consider concave aberration criteria, some of our results still hold. 

Consider any concave function / : R d — ► R. Now, given a model L, define 
its aberration by 

A(f,L) := f lj2 a i,-iJ2 ad ) ■ 

\aeL a£L / 
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The linear aberration of Definition [T] is the special case where / is the fol- 
lowing linear (hence concave) function, 



./ 



x 



(xx, ...,x d ) 



1 d 

n ^ 



WiXi 



Since we only appealed to convexity, Theorem [3] is valid when we replace 
A(w,L) by the more general form A(f,L). That is to say, the set of lower 
vertices of the state polytope (corresponding to models in the algebraic fan) 
contains the solution to minimising any concave aberration function. This 
can be understood as minimisation over a matroid, which was studied further 



in 



Berstein et al.l (120081 ). A further development is to consider aberration 



A(w, S(w)) with respect to other distributions rather than the uniform. 



5.2 Connection with aberration of Wu and others 

In the statistical literature, the word aberration has been used to refer to 



prope rt ies of the generato r s for f ra ctional factor i al des igns, see I Chen and Hedayat 



(11998|); iFries and Hunted fll980f ): IWu and Wul ((20021). A topic of future re- 



search is to link minimal aberration of Definition [T] with the traditional mea- 
sure based on generators for a fractional factorial design. 

We conjecture that among the class of orthogonal fractions of 2 d designs 
there is some kind of correspondence between the minimal linear aberration 
of this paper and minimum generator aberration of Wu and others. If we 
select non-orthogonal fractions, the situation is more complex, as the next 
example shows. 

Example 10 Let F 3 be the non-orthogonal fraction with size n = 32 of a 2 7 
design given in Table [2] of Appendix B. We also consider the designs Fi and 
F 2 of Examples [8] and [9j The three designs have the same size, but the design 
F 3 cannot be compared with Fi or F 2 in traditional terms as it is not even 
orthogonal. However, we can compare the designs based in the distribution 
of degrees in their algebraic fans. 

An interpolation as presented in Appendix B suggests that the mini- 
mum degree of models identified by a generic design with n = 32, d = 7 is 
53.5 ~ 54. This number is a lower bound for the total degree of models 
identified by designs F 1 ,F 2 and F 3 . In other words, the set of total degrees 
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for models in algebraic fan of Fi, F 2 and F 3 is lower bounded by 54, e.g. 
54 < min({£ti a L : L G £«(*!)}) for i = 1, 2, 3. 
Initial results show that 



i) the size of £ a (i<3) is much longer (it has around 6 x 10 5 models) than 
that for designs F% and F 2 , see Tabled] in Appendix B; 

ii) the algebraic fans of Fi and F 2 are not contained in the algebraic fan 
of F 3 , and 

iii) the design F 3 identifies model of lower degree than Fi or F 2 (indeed of 
total degree 58), and the bound 54 is verified. 

It is clear that F3 has smaller minimal linear aberration than F\ and F 2 , see 
Figure [121 We also note that the histogram for F3 presents more symmetry 
than that for Fi and F 2 . 
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Appendix A: Proof of Theorem H 

Proof. The proof is basically the minimisation over the first orthant of 
Ylt=i w i x i subject to the constraint n^il 3 -* + a) — b. The problem is solved 
by a change of coordinates to i- = Xi + a for i = 1, . . . , d. We minimise 
Ylt=i w i x 'i subject to Yli=i x[ — b = . Using standard optimization tools, 
we form the Lagrange multiplier 
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and then solve the system of equations VL(x', A) = 0, — = 0. The 
solution vector *, ) where 

T-rd 1/d 

xf = b l ' Aii = lWl . 

Wi 

The convexity of the functions Yli=i w % x i an d EliLi x i = b over the nrs t 
orthant guarantees that x* is indeed the minimum. The aberration for this 
minimal point is 

d 

^WiX* = db 1/d g(w). 

i=l 

Finally we note that x* = x* — a and compute the aberration using x*, 
achieving the approximate aberration A of Equation ([7]). ■ 

We remark that for a fixed w, x* serves as an approximation to the cen- 
troid of the corresponding corner cut model and therefore A is an approxima- 
tion to A*(w, n). Although the approximate aberration A does not depend on 
the actual corner cut identified by L, the minimal aberration A*(w,n) does 
depend on it. If L is the corner cut directed by w, the practical validity of 
the approximate aberration A relies on x* being close enough to ^ ^2 aeL 
This closeness depends ultimately on a, b. See Appendix 15.21 for a proposal 
to compute a, b. 



Appendix B: Computing values a, b for the ap- 
proximate corner cut polytope 

In Section 13.41 we proposed the continuous function of Equation ([6]) to ap- 
proximate the corner cut polytope (which is piecewise linear surface). In this 
section we discuss on the selection of the values a, b so that the approxima- 
tion is good enough. In general, the values a, b will depend on the number 
of dimensions d and number of points in the design n. However, for fixed d, 
the approximation will be coarse for small values of n. 

For our approximation we use the follo wing properti es of the corner cut 



polyto pe, which have been studied as well in lMullerl (120031 ) and lOnn and Sturmfels 



fll999[ ). 



Lemma 10 The corner cut polytope satisfies the following properties. 
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i) The intersection of the corner cut polytope with the axes occurs at the 
point (™) . 



ii) When for k > 1, the sample size n satisfies 

(9) 



n 



k + d-1 
d 



then the corner cut polytope is pointed. 
Proof. 

i) The intersection is the the sum of exponents for any marginal model 
of the form {1, Xi, x?, . . . , a;™ -1 }. Therefore the intersection must occur 

at Yh=o * = (2) ■ 

ii) The corner cut polytope is pointed when the sample size is the same 
as the size of a model of total degree k — 1, that is, there are 
terms of degree j in the model where j — 0, . . . , k — 1. Therefore the 
sample size must be n = ^tj = (^V)- 



Remark 11 When Equation (GJ) is satisfied, the tip of the pointed corner 
cut polytope has coordinates ai = ({ k ^~[ 1 ), • • • ' (^d+i" 1 )) ■ 

We propose to force Equation (jSj) to satisfy the condition of Item 1 in 
Lemma [TU1 and pass through the tip point ai for the model of total degree 
k — 1. To summarize, when sample size satisfies Equation dHJ) then a, b must 
satisfy the following equations: 



n 



and b = (c + a 



where c = ^ ( fc ^ x 1 ) is the scaled tip of the corner cut poytope. When design 
size, n, is not of the form n = ( fe+ ^ _1 ) for some k > 1, we propose to 
interpolate the value for c, the scaled tip of the polytope, that is to solve 
Equation ([9]) for k and interpolate the corresponding tip with ^{ k ^ l+ ^ 1 )- 

For two dimensions (d = 2) by interpolation and solving the two condi- 
tions above we obtain the following formulae for a, b in terms of n: 



5 - 3V1 + 8n + An , fn-l 

a = — -, b = a \ — h a 

3(3 - 2-n/I + 8n + 3n) 
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See Figure d3] for a depiction of the corner cut polytope and the approximate 
curve for d = 2, n = 7. This interpolation is difficult for d > 2 and we have 
to rely on approximations. The following formula? are rough approximations 
for a, b obtained by truncation of the binomial expansions 




1 2 ' ^3 



Figure 13: Minimal aberration using the corner cut polytope. The corner cut 
polytope is the piecewise linear solid curve, while the approximation is the 
dashed curve. The minimal aberration is the projection over the direction of 
w of the vertex (using dotted line), and an approximate value uses Equation 
(El) (dashed line). 
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Total 
degree 


AF F± 


AF F2 


AF F3 


RF Fi 


RF F2 


RF F3 


58 


- 


- 


2290 


- 


- 


0.84 


59 


- 


- 


5437 


- 


- 


1.99 


60 


- 


- 


15036 


- 


- 


5.51 


(il 


- 


8 


34574 


- 


0.47 


12.66 


62 




52 


55025 




3.04 


20.15 


63 




108 


57848 




6.32 


21.18 


64 




124 


47851 




7.26 


17.52 


65 




220 


28511 




12.88 


10.44 


66 




268 


13928 




15.7 


5.1 


07 




204 


6837 




11.94 


2.5 


68 


72 


340 


3378 


54.14 


19.91 


1.24 


69 




60 


1596 




3.51 


0.58 


70 




136 


567 




7.96 


0.21 


71 




8 


140 




0.47 


0.05 


72 


48 


144 


33 


36.09 


8.43 


0.01 


73 






12 






0.00 


74 




20 


5 




1.17 


0.00 


8(1 


12 


16 




9.02 


0.94 


0.00 


83 






1 






0.00 


85 






1 






0.00 


119 


1 






0.75 






Total 


133 


1708 


273071 


100.00 


100.00 


100.00 



Table 1: Absolute (AF) and relative (RF) frequencies of total degrees for 
models identified by fractions F\ and F% of Example [9] and F3 of Example 
[TUl The symbol - represents zero. 
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